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I Abstract. Let a„ denote the number of ways in which a chess 

Qh' rook can move from a corner cell to the opposite corner cell of 

an n X n X n three-dimensional chessboard, assuming that the 
piece moves closer to the goal cell at each step. We describe the 
I computer-driven discovery and proof of the fact that the generating 

series G{x) = J2n>o ^n^" admits the following explicit expression 
in terms of a Gaussian hypergeometric function: 
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^ ■ 1. Introduction 

in ; 

^ ■ In this article, we solve a problem of enumerative combinatorics ad- 

dressed and left open in [16]. The initial question is formulated in terms 
^ \ of paths on an infinite three-dimensional chessboard. The 3D chess- 

board being identified with N^, a 3D rook is a piece which is allowed 
to move parallelly to one of the three axes. The general objective is 
to count paths (i.e., finite sequences of moves) of a 3D rook on the 3D 
^ \ chessboard. Following [16], we further restrict to 3D rook paths start- 

ing from the cell (0, 0, 0), whose steps are positive integer multiples of 
(1, 0, 0), (0, 1, 0), or (0, 0, 1). In other words, we assume that the piece 
moves closer to the goal cell at each step. For example, one such rook 
path is 

(0,0,0) ^ (5,0,0) ^ (5,0,6) ^ (10,0,6) ^ (10,2,6). 

Let Tj ,,• denote the number of rook paths from (0, 0, 0) to (z, j, /c) G 
N^, and let a„ = rn^n.n be the number of diagonal rook paths from 
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(0,0,0) to {n,n,n). The sequence (a^) is called the 3D diagonal rook 
sequence; its first terms are 

1, 6, 222, 9918, 486924, 25267236, 1359631776, 

75059524392, 4223303759148, . . . 

This is Sequence A144045 in Sloane's on-line encyclopedia of integer 
sequences at http://oeis.org/A144045. It was conjectured in [16] 
that the sequence (a„) satisfies the fourth-order recurrence relation 

(1) 2n'^{n - l)a„ - {n - l)(121ra2 - 91n - 6)a„_i 

- (n - 2) (475^2 - 2512n + 2829)a„_2 + 18(n - 3) (97n^ - 519n + 702)a„_3 

- 1152(n-3)(n-4)2a„_4 = 0, forn > 4. 

The aim of the present article is not only to give a computer-driven 
proof of this conjecture, but also to show that the generating function 
G{x) = Z]n>o'^n^" ^ admits the following exphcit expression: 
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where the Gaussian hypergeometric series 2-^1 with parameters a,b,c E 
C, — c ^ N, is defined by 

a b 
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using the notation (a)„ for the Pochhammer symbol (a)„ = a{a + 
l)---{a + n-l). 

The core of this article is Section 2, in which we complete the enu- 
merative study of 3D rooks, proving in particular recurrence (1) and 
the explicit form (2). Section 3 then presents further comments, variant 
proofs and calculations, and discussions of a possible extension. 

In more detail. Section 2 proceeds as follows. The question of the 
enumeration of a constrained combinatorial walk is re-expressed in Sec- 
tion 2.1 in terms of a diagonal of a rational series. Then, this diagonal 
is encoded as a residue in Section 2.2, leading to the problem of finding 
a special ODE it satisfies. In the literature, an effective approach to 
finding this special ODE was developed by Lipshitz. In Section 2.3, 
we show that, while the approach applies in theory, it leads to unfea- 
sible calculations in practice. A faster approach, Zeilberger's creative 
telescoping and its evolution by Chyzak, is summarised and applied 
in Section 2.4. The obtained ODE allows us to prove the conjectured 



EXPLICIT FORMULA FOR DIAGONAL 3D ROOK PATHS 



3 



recurrence (1) and even a lower-order one, (15), in Section 2.5. We 
then derive the 2-^1 exphcit form in Section 2.6. 

When possible, we have attempted to describe the mathematical 
phenomena, ideas, and techniques on a general level in Sections 2.1- 
2.6. When the complete proofs of the results we state require a more 
involved calculation, and also when they are based on more involved 
notions or on the fine tuning of algorithms, we have postponed some 
of the technicalities to Appendix A. 

2. Main steps of the solution 

2.1. Prom combinatorics to algebra. The combinatorial problem 
translates into a question on generating series, best expressed in terms 
of diagonals: given (f){s,t,u) = J2i,j.k>oCi,j,ks''Pu'' in Q[[s, t,^]], its diag- 
onal is defined as the univariate power series Diag(0) = J2n>o Cn,n,nX"' 
in QWx]]. 



Question. Define the trivariate rational function f{s, t, u) in Q(s, t, u) 



as 



n>l n>l n>l 



s){l-t){l-u) 



l-2{s + t + u) + ?>{st + tu + us)- istu ' 

which is seen to be at the same time in Q[[s, t,n]]. What can be said 
about its diagonal Diag(/) = J2n>o '^nX^ ? 

This translation is based on the following straightforward extension 
of Proposition 6.3.7 in [31], which we shall use with the choice (3 = 
{ (n, 0, 0) : n > 1 } U { (0, n, 0) : n > 1 } U { (0, 0, n) : n > 1 }. 

Lemma 1. Given a finite set of allowed steps © C \ {(0, 0, 0)}, let 
Ne{m,n,p) denote the number of paths from (0,0,0) to {m,n,p), with 
path steps taken from &. Then, the generating series 

Ge(s,t,u)= N6{m,n,p)s"'euP 

m,n,p>0 

is rational, and has the form 

G6{s,t,u) = i 

The same result holds for infinite step sets & provided the power series 
J2{i,j,k)ee s^t^u'' is rational. 
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2.2. The key equation. A first, theoretical answer to our question 
comes from the main result in Lipshitz's article [24]: the diagonal 
Diag(/) is a D- finite series, i.e., it satisfies a linear differential equation 
in dx = d/dx with polynomial coefficients in x. Therefore, the sequence 
{an) of its coefficients satisfies a linear recurrence with polynomial co- 
efficients in n. 

Moreover, let F be the trivariate rational function 



F=--f{s,t/s,x/t) 
st 



[s - t){s - l){t - x) 



st{~st + 2sH + 2t2 + 2xs - 3st2 - 'ixt - Sxs^ + Axst) ' 
Then, Remark 3 in [24] shows that 
Lemma 2. // the equation 

admits a solution {P,S,T) where P{x,dx) is a non-zero linear differ- 
ential operator with polynomial coefficients in Q[x], and S and T are 
two rational functions in Q{x,s,t), then P{x,dx) annihilates Diag(/). 

Lemma 2 is the heart of the result in [24], and the basis of all methods 
for effectively computing differential equations satisfied by diagonals. 
For completeness, we give here a self-contained proof, which closely 
follows [24]. 

Proof. Although F is not a formal (Laurent) power series in x, s, t, it 
does come from a formal power series after a change of variables, and 
thus it can be seen as an element of the differential Q[x, s, t]-module of 
all formal sums of the form 



B = J2 Ca,l3,'yS"t'^X 



i + ^>-M,/3+7>- 



for some (yU, u) G depending on B. As a consequence, it makes sense 
to speak of the series expansion of F and of P{F), but in later calcu- 
lations, these series may be multiplied solely by elements of Q[x, s,t]. 

The diagonal Diag(/) is, by construction, equal to the coefficient 
of s~^t~^ in F. On the other hand, equation (3) implies that the 
coefficient of s~^t~^ in P{F) is zero. The conclusion then follows from 
the sequence of equalities 

P(Diag(/)) = P([."V^]F) = [s~H-']P{F) = 0, 
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the second equahty being a consequence of the fact that P is free of dg 
and dt- □ 

The question remaining is how to determine a non-trivial triple 
{P,S,T) solution of (3). 

2.3. Lipshitz's approach. Lipshitz [24, Lemma 3] uses a counting 
argument, similar to Fasenmyer's [17, 39], to show that there exists a 
non-zero linear differential operator L[x, d^, dg, dt), free of s and t, that 
annihilates F. Then, writing this operator under the form 

L(x, d^, ds, dt) = P(x, dr,)+ 

(higher-order terms in dg and dt, free of s and t), 

it follows from Lemma 2 that P annihilates Diag(/). 

Let us sketch Lipshitz's argument on our concrete example; as we 
shall see, while this proof is algorithmic and ultimately reduces the 
existence of P to a linear-algebra argument, it produces unreasonably 
large systems. The key observation is that, by Leibniz's rule, for any 
nonnegative integer A^, the (^^^^^ rational functions 

x'did% (F), <i+j + k + i < N, 

belong to the Q- vector space of dimension at most 18(A^+ 1)'^ spanned 
by the elements 

(4) < z < 2iV + 1, < J <3N + 2, < k < 3N + 2, 

where q is the denominator st{—st + 2s^t + 2t'^ + 2xs — 3st'^ — 3xt — 3xs'^ + 
Axst) of F. Thus, if is an integer such that (^+^) > 18(iV + 1)^, 
then there exists a non-zero linear differential operator L{x,dx,ds,dt) 
of total degree at most in x, d^, dg, and dt, such that L{F) = 0. 

Now, it appears that = 425 is the smallest integer satisfying 
the inequality > 18(A^ + 1)^; therefore, exphcitly finding the 

operator L would require solving a homogeneous linear system with 
some 1,391,641,251 unknowns and some 1,391,557,968 equations! 

It is possible to refine the previous argument in various ways, and 
thus to reduce the search of L to a linear algebra problem of smaller 
size. Such methods will be discussed and compared thoroughly in a 
later work. In our case, the best reduction in size that we could achieve 
uses a variation of the main idea in [9, §3.2]. 

First, as before, Leibniz's rule shows that the un '■= {^^^^ rational 
functions 

did^d'tiF), 0<j + k + i<N, 
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belong to the Q(x)-vector space spanned by the elements 

(5) < J < 3iV + 2, 0<k<3N + 2, 

whose dimension djsf is at most 9{N + 1)^. 

A further idea is to exploit the sparse feature of the monomial sup- 
ports of the numerator and denominator of F. By studying the be- 
haviour of these supports under differentiations with respect to x, 
s, and t, one shows that the terms (5) will not be involved for the 
^2N+2^^ + (^+2^ pairs (j, k) e N2 in the set 

{{j,k):j + k<2N} U 

{ (j, k):2N + 2<j,k<3N + 2, j + k>5N + 4}. 

This implies that dN is bounded above by cn '■= |(13A^-|- 14)(A^+ 1), 
and shows that the search of L boils down to a homogeneous linear 
system over Q{x), with unknowns and ejy equations. The smallest 
value of N such that mat > cat is iV = 36, leading to the task of 
computing the (right) kernel of a matrix of size 8917 x 9139, with 
entries in Q[x] of degree at most 37. 

Despite the spectacular reduction in dimension compared to Lip- 
shitz's original argument, these sizes are still too high to be dealt with 
in practice. Of course, the value cn is only an upper bound on dN, 
leaving room for a magical rank drop occurring before the predicted 
value = 36. That being said, our implementation of the method 
constructs the linear system for increasing values of N, but does not 
find any non-trivial kernel for N < 18, and is not able to finish the 
computation for = 19 before consuming all the available memory.^ 

The conclusion is that Lipshitz's approach and its improvements 
are not sufficient to obtain effectively an equation for the diagonal 
of /(x, 

2.4. The creative-telescoping approach. An alternative approach 
to solving (3) was popularised among the computer-algebra commu- 
nity by Zeilberger under the name of creative telescoping, and it is 
the one that we shall use to effectively construct an equation for our 
diagonal. This method, first introduced for hypergeometric summa- 
tion [40, 42] and hyperexponential integration [3], was later generalised 
by Chyzak [11] to more general integrands involving a larger class 

^Our implementation of an optimisation of Lipshitz's algorithm, available at 
http : //algo . inria.fr/ chyzak/Rooks/OptimisedLipshitz .mpl, is nevertheless 
able to predict that, with high probability, a non-trivial kernel does exist for — 19. 
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of special functions called holonomic. Zeilberger's definition of this 
class [41, Section 2.2.4] has its roots in D-module theory [5, 6] and is 
not needed to understand the present paper. See also [12] for recent 
extensions of creative telescoping to more general special functions. 

Prior to summarising the method, we show its behaviour and its re- 
sult on our concrete case, postponing the explicit details of intermediate 
steps to Appendix A. 

2.4.1. Obtaining the result. In our case, Chyzak's algorithm, as imple- 
mented by the package Mgfun^ for Maple, produces a solution (P, S, T) 
of (3), with P = P2dx, where 

(6) P2 = x{x - l)(64x - l)(3x - 2)(6x + l)dl 

+ (4608x^ - 6372x^ + 813x2 + 514x - 4)9^ 

+ 4(576x3 - SOlx^ - 108x + 74). 

Let q = st ■ qi with qi irreducible in Q[x,s,t] of degree (1,2,2) in 
(x, s, t). Then, S and T are of the very special forms: 

c is~t)-U ^ {s-t)-V 

2st ■ qf ■ disct(gi) ' 2s'^ ■ qf ■ disCf(gi)^ ' 

where U and V are (irreducible) polynomials in Q[x, s,t] of respective 
degrees (5,8,3) and (8,14,5) in (x, s,t). Here, disCf(gi) denotes the 
discriminant of qi with respect to t. Observe the very moderate degrees 
in comparison to those pessimistically predicted by Lipschitz's method. 

Maple commands to produce P, S, and T are provided online at 
http : //algo . inria.fr/chyzak/Rooks/CreativeTelescoping.mpl; 
they take less than 10 seconds on an average laptop. The rational 
functions S and T are not displayed here, but can be found online at 
http : //algo . inria. f r/chyzak/Rooks/Certif icates .mpl. 

2.4.2. Sketch of the algorithm. We proceed to sketch the algorithm 
of [11] for the creative telescoping of multiple "integrals." By this, we 
actually mean in this article answering the algebraic question of how 
to solve the key equation, whether in the form (7) below for "simple 
integrals," or in the form (3) for "double integrals." 

Presenting Chyzak's algorithm directly in its generality would prove 
too arid, so we go instead progressively from the simpler case of rational 
simple integrals to the integration of more general functions, then to 
double integrals. In effect, this section is mostly bibliographic, with the 

^version 4.0, available as part of Algolib 13.0 from 
http://algo.inria.fr/libraries/. Initially developed by F. Chyzak. Specific 
version used here reimplemented by L. Pech. 
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exception that, in the case of multiple integrals, we provide an extension 
of the algorithm in [11] beyond the so-called "natural boundaries." 

For a single integral of a bivariate rational function F{u,v), the 
creative-telescoping method finds a linear differential operator P{u, dy) 
and a rational function S in Q{u,v) such that (cf. equation (3)) 

(7) Pin - f . 

To this end, the algorithm in [3] writes P and S in undetermined form 

(8) P = r]r{u)dl + ■ ■ ■ + r]o{u) and S = (j){u,v)F 

for some tentative order r and unknown rational functions r)i from Q{u) 
and from Q(m, v). To get 0, the original algorithm of [3] exploits the 
specific nature of F (hyperexponential). We prefer describing here 
the less optimal variant of [11], as it will extend to the more general 
class of inputs we need: (7) is rewritten into an auxiliary ordinary 
linear differential equation on 0, before one uses known algorithms for 
solving a linear differential equation with polynomial coefficients for its 
rational solutions. More precisely, such algorithms, like the one in [1], 
have to be refined to solve for the rji as well; this is easily taken into 
account just by having more unknowns in the crucial step of linear 
algebra over Q{u) in this rational solving. 

For reasons to be given shortly, the double integration of trivariate 
rational functions is based on integration of more general functions F, 
leading to a more general form for S in (7). Specifically, the case dealt 
with in [11] is that of functions F for which there exists a finite maximal 
set of derivatives d'^dl{F) that are linearly independent over Q{u,v). 
Such functions are known as differ entiahly finite in the literature, or 
in short D-finite [30, 25]. However, note that the notions of holonomic 
and differentiably finite functions, which look technically very different, 
cover in fact exactly the same sets of functions. This fact is explained 
in elementary terms in [32, Appendix]. Still, the more elaborate holo- 
nomic point of view is needed in proofs that require a subtler counting 
argument. 

Once such a set of pairs (a, h) is fixed for a given D-finite function F, 
equation (7) is changed by setting S to the undetermined form 

(9) S=Y.<PaAn.v)d:dl{F) 

{a,b) 

for unknown rational functions 0a,6 from Q(m, t"). The only change 
then is that the auxiliary equation in becomes an auxiliary ordinary 
linear differential system in the (f)a,bS. Solving is performed either by 
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uncoupling and using the algorithm mentioned previously several times, 
or by a direct approach as in [4]. 

The case of double integrals leads to generalising (7) into (3), but 
the solving does not generalise so smoothly: attempting to write S = 
(j)i{s,t,x)F and T = (j)2{s,t,x)F for undetermined rational functions 
01 and 02 from Q(s,t,x) yields a linear partial differential equation 
relating these functions with d(j)i/ds and 902 To the best of our 
knowledge, although this overdetermined equation has a very specific 
form, no algorithm is available to solve it for its rational solutions. 

Therefore, instead of a direct approach, Chyzak developed in [11] 
a cascading approach which we now summarise. Noting that the de- 
pendency of P on a single derivation du in (7) is inessential, the same 
approach is possible for the creative telescoping with respect to the 
(single) variable i; of a trivariate rational function F from Q{ui, U2, v). 
Indeed, setting P to the undetermined form 

(10) P= i: m,A^i.n2)dix^ 

0<i+j<r 

for some tentative total order r and unknown rational functions rjij 
from Q{ui,U2) and performing the same solving as previously, now 
relying on linear algebra over Q{ui, U2), leads to a basis of P*^°^'s of total 
order at most r for which there exists a rational function 0^°)(mi, U2, v) 
satisfying 

p(°)(p) = ^ f0(")F) . 

ov ^ ' 

The theory (as developed, e.g., by Zeilberger in [41]) guarantees that 
the function F is D- finite and that the set of P^^^'s obtained for suf- 
ficiently large r can be used to determine the finite set needed as an 
input to the algorithm of [11] and in (9). 

Finally, a double integration algorithm is obtained by continuing the 
approach used for natural boundaries in [11] (Stages A and B below) 
by a suitable recombination of the outputs (Stage C below). The re- 
sulting treatment of multiple integrals over non-natural boundaries is 
an extension over [11], and the corresponding algorithm is as follows: 

• Stage A: First iteration of creative telescoping. Using the uni- 
variate algorithm for trivariate rational functions and variables 
(s, a;,t) in place of (ui,U2,w) delivers identities 

(11) P(°)(., x, 9„ d,){F) = ^ (0(")(., t, x)f) . 

• Stage B: Second iteration of creative telescoping. Considering a 
function F of (s, x) that is annihilated by all P*^"-* and using the 
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univariate algorithm for general functions and variables (s,x) 
in place of {u, v) delivers an identity 

(12) P{x,d^){F) = ^^[Q{s,x,ds,d,){F)). 

• Stage C: Recombination. By the theory of linear-differential- 
operators ideals,'^ the calculations of the algorithm can be in- 
terpreted as a proof of existence of operators L^°'\s,x,ds,dx) 
satisfying 

(13) P(x, d.,)-dsQ{s, X, ds, dx) = J2 L^"'^ (s, X, ds, 9,.)P("ns, a;, 9„ d^). 

a 

These L^^^ can effectively be obtained either by following the 
calculations step by step or (less efficiently) by a postprocessing 
(non-commutative multivariate division). Hence, defining 

(14) 

S = Qis,x,ds,d,)iF) and T = ^ L(")(s, x, 9„ 9,) (s, t, x)f) 

a 

leads to a solution (P, S,T) of (3). 

Note that this two-stage process inherently introduces a dissymmetry 
in the treatment of the variables s and t: the output from the first 
iteration tends to be larger than its input; in turn, the output from the 
second is larger than the output from the first. As a consequence, the 
order we deal with the variables may have an impact on the running 
time. 

2.5. Proving the recurrence. Translating the differential equation 
P(Diag(/)) = into a recurrence satisfied by the coefficients a„ of 
Diag(/) yields the fourth-order recurrence relation (1). 

It is possible to prove that the sequence (a„) satisfies the even shorter 
recurrence 

(15) 

2 (n - 1 ) (35n - 52)n2a„ - (n - 1) (4655n^ - 1 1 78 In^ + 8494n - 1 776) a„_i 
+ {n- 2)(11305n3 - 41856^2 + 46487r2 - 13128)a„_2 

- 192(n - 3)^(35/2 - 17)(n - 2)a„„3 = 0, for n>3. 

Recurrence (15) can be guessed (e.g., by linear algebra) using, say, 
the first 25 terms of the sequence (a„). Once it is conjectured, proving 

"'This observation was made during a discussion between F. Chyzak, M. Kauers, 
and Ch. Koutschan on July 17, 2008, in the "UFO" room at RISC (Linz, Austria). 
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its correctness is easy: if denotes the sequence defined by tlie left- 
hand side of (15), then recurrence (1) imphes that {bn) satisfies the 
recurrence 

bn + 66„_i = {35n - 52) (2{n - l)n'^an -{n- 1) (121^2 - 91n - 6)a„_i 

- (n - 2) (475^2 - 2512n + 2829)a„_2 + 18(n - 3) (97^^ - 519ra + 702)a„_3 

-n52{n-3){n-4:fan-?j=0, for n > 3, 

which, combined with the initial condition b^ = —54864 ■ 6 + 2 ■ 43362 ■ 
222 - 2 ■ 2 ■ 53 ■ 32 ■ 9918 = 0, yields 6„ = for all n>0. 

Maple code to obtain the first terms of the sequence and derive the re- 
currences can be obtained at http : //algo . inria.fr/chyzak/Rooks/Recurrence .mpl. 

2.6. Solution in explicit form. Operator P2 from equation (6) has 
a power series solution with integer coefficients, namely c?2:(Diag(/)). 
By a speculation of Dwork's [15], it appears plausible that P2 has a 
2-Fi-type solution. We shall briefly sketch how, for a second-order op- 
erator L, one might find such solutions. A similar problem (for Bessel 
instead of 2F1) was solved in [14, 36], and many of the details discussed 
in those papers are needed for the 2-^1-case as well (in particular the 
first part of Section 3 in [14]). The result of our computation is: 

'"""-'™- (l-4.)(l-64.) -^'C'"2''" 



27a;(2-3x)^ 
(1 -4x)3 



2.6.1. Formal asymptotic- series solutions of linear differential opera- 
tors. We recall the necessary facts from the classical theory of asympto- 
tic-series solutions of linear differential operators [22, 38], and specialise 
them here to our need: the order 2. 

With respect to a given linear differential operator L of order r with 
polynomial coefficients, each point p from C U {00} falls in one of two 
categories: 

• it is called an ordinary point of L if the coefficient of the highest 
derivative in the equation does not vanish at p, in which case 
a solution / is fully determined by the initial conditions f{p), 
. . . , f^^~^\p) (Cauchy theorem); 

• it is called a singular point of L otherwise. 

Among singular points, the more tractable situation is that of a reg- 
ular singular point (of L). At such a point, the growth of solutions is 
bounded (in any small sector) by an algebraic function. In the neigh- 
bourhood of a regular singular point p, any solution can be written as a 
linear combination of products of a (possibly multivalued) term of the 
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form {x — pY with e G C, logarithms ln{x—p), (respectively {1/xY and 
ln(l/x) ii p = oo) and analytic (univalued) functions. The constants e 
are called the exponents of L a,t x = p. They are the zeroes of the 
so-called indicial equation of L at x = p, an algebraic equation that is 
defined for any L and p, and can be obtained algorithmically. 

Another simple situation occurs at a singular point p of L when 
there exist a function u (possibly singular at p) and functions yi, . . . , 
Ur that are analytic in a neighbourhood of p, such that uyi, . . . , uyr 
constitute a basis of the solution space of L. Then, we say that p is a 
removable singularity of L. Otherwise, we say that it is a non-removable 
singularity. 

Assume from now on that L is irreducible and of order 2. Then, the 
indicial equation a.t x = p has degree at most 2; in other words, there 
are at most 2 exponents ai x = p (counting with multiplicity). We 
shall also assume that L is regular singular ai x = p, which implies 
that the degree of the indicial equation is equal to 2 and that there are 
two exponents, Ci and 62 (counting with multiplicity). The exponent 
difference d at x = p is defined up to sign as ±(ei — 62). 

The numbers ei, 62, and d describe a lot of the asymptotic behaviour 
of the solutions of L when x tends to p. Let us be more specific about 
this. For notational convenience, consider first the case p = 0. Then 
e is an exponent at x = if and only if L has a formal solution u for 
which u/x'^ G R\xR, where R := C[[x]][ln(x)]. If ln(x) actually occurs 
in u/x"^ then the point p is called logarithmic. If ln(x) does not occur, 
then our assumption that L is regular singular can be used to prove 
that u/x^ is analytic at x = 0, so that x^ describes the local behaviour 
of M at X = 0. Similar statements hold for p ^ 0, after replacing x 
by X — p if p is finite, or x by 1/x if p = 00. The following lemma 
summarises the link between the exponent difference of a second-order 
operator and the asymptotic behaviour of its solutions. 

Lemma 3. Let L be a linear difference operator of order 2 with expo- 
nent difference d at some point j> G C U {00}. Assume p to be either 
an ordinary or a regular singular point of L. Then: 

(1) If p is ordinary, then d = 1. 

(2) If p is a removable singularity, then (i G Z \ {0}. 

(3) If p is logarithmic, then d ^'L. 

(4) If d = 0, then p is logarithmic. 

(5) If d G Z \ {0}, then p is either logarithmic or a removable 
singularity. 



For proofs, see [14, Theorem 1], or see [13] for more details. 
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2.6.2. Transporting the singularities of the Gauss hypergeometric equa- 
tion. Let L^'^ denote the differential operator of the Gauss hypergeo- 
metric equation, 

x[l — x)—-Tr + c — (a + + Ijx — abw = 0. 

ax'^ ^ ^ ax 

Of its two hnearly independent classical solutions, we shall only need 

one: 2-^i(«5 b; c; x). It has 3 singular points: 0, 1, and oo. The respective 

exponent differences are (eo,ei,eoo) = (±(1 — c), ±(c — a — 6), ±(a — 6)). 

Let / be a rational function, and denote by the operator obtained 

from L^'^ by the change of variables x = f, so that it has 2-^1 (cf? b; c; f) 

as a solution. 

The following lemma explains how the exponent differences of the 
Gauss differential equation get transported by the composition under /. 
The proof is the same as in [14, Theorem 2(i)] (more details are given 
in [13]). 

Lemma 4. Let d he the exponent difference of L'^'.^'f at x = q. Then: 

(1) If q is a root of f{x) with multiplicity m, then d = m ■ cq. 

(2) If q is a root of f{x) — 1 with multiplicity m, then d = m ■ Ci. 

(3) If q is a pole of f{x) of order m, then d = m ■ Coo- 

Recall that L is irreducible and of order 2. Let y = 2-P'i(a, 6; c; /) 
where / is some rational function. By solving L in terms of y, we mean 
finding three rational functions r, tq, and Vi such that 

(17) exp (^J ■ {roy + riy') 

is a solution of L (see [14, Definition 2, parts (ii) , (Hi)]) . A neces- 
sary condition for such a solution to exist is that L has the same set 
of non-removable singularities as -Z^c-p moreover, that the expo- 
nent differences match modulo Z (see [14, Lemma 4]). The expression 
exp(/ r) is usually a radical function, i.e., a function whose rth power 
is rational for some integer r. 

Consider the operator P2 in (6). Its non-removable singularities are 
0, 1, 1/64, 2/3, and 00, and all these singular points are logarithmic; 
on the other hand, the root x = —1/6 of the leading coefficient of P2 
is a removable singularity. 

Suppose that we choose a,b,c so that (cq, 61,600) = (0,1/3,0) (this 
choice will be explained in Section 2.6.3), and want to solve P2 in terms 
of y := 2F'i{a,b; c; f). By the discussion above, if p is a logarithmic 
point for P2 then it must be a logarithmic point for and vice versa. 
Since x = and x = 00 are the logarithmic points of L^'^ one sees 
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that the union of the set of roots and the set of poles of / must be 
exactly {0, 1, 1/64, 2/3, oo}. For each possible degree of /, this leaves 
finitely many candidate integer values for no, ni, ni/64, and 7^2/3, such 
that / is of the form c(x - 0)"o(x - l)'''{x - 1/64) '"1/64 (2; _ 2/3)''^^ 
for an unknown constant c. For an exhaustive search, we therefore try 
successive potential degrees 2, 3, 4, . . . for /. 

The exponent difference of at x = g is m ■ ei if g is a root of / — 1 
with multiplicity m. So m ■ ei must match (modulo Z) an exponent 
difference of P2. Since P2 has only integer exponent differences, and 
Ci = 1/3, it follows that 3 divides m and so the numerator of / — I must 
be a cube. Thus, we try various tt-q, ni, . . . G Z, and for each, we check 
if there exists c G C for which the numerator of c{x — 0)""(x — l)"i(x — 
l/64)"V64(a;_2/3)"2/3_i is a cube. A computer implementation quickly 
finds 

no = 1, ni = -2, 122/3 = 1, ?^i/64 = -1, c = -81/64, 
and we get 

-81x{x - 2/3) 
~ 64(a;- l)2(a:- 1/64)' 

Next, we compute a,b,c for which the exponent differences eo,ei,eoo 
are 0,1/3,0 and find L^'''. Applying a change of variables produces 
L'^j. Then a call to the program equiv [34] provides a map of the form 

y (— ^ exp(/ r)(ro?/ + riy') from the solutions of L".j to the solutions of 
P2. This way we find a solution of P2 in the form (17). 

2.6.3. Simplifying the obtained solution. The question now is how to 
make this solution more compact. First of all, the numerator of / — 1 is 
a cube. Interchanging ei and 600 makes the denominator of / a cube. 
That makes the expression for / slightly smaller. Now, composition 
with x/{x — 1) interchanges 1, 00 while keeping fixed, so we replace 
/ by // (/ — 1) and obtain a new /, 

_ 27a;(2 - 3a;) 
^ ~ (l-4x)3 ' 

which corresponds to the new (eo, 61,600) = (0,0,1/3). Next, we try 
to change 60, 61, 600 by some integers in such a way that the exponent 
differences of not only match those of P2 modulo Z, but are also as 
close as possible, because this tends to reduce the size of ro, ri. Surpris- 
ingly, one can even match the exponent differences exactly! That allows 
ri to become 0. We get an exact match for (60,61,600) = (1,1,1/3), 
corresponding to a = 1/3, 6 = 2/3, c = 2, and end up with the following 
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solution 



6 -,F^^^ 



;i-4x)(l-64x) 



27x{2-3xy 
(1 - 4x)3 



This expression is not unique because we can interchange (eo,ei,eoo), 
updating / accordingly. We can also change eo,ei,eoo by integers (as 
long as those integers add up to an even number) but doing that would 
make ri non-zero again. Moreover, there are other exponent differences 
that would have solved P2 as well. 

The set of (eo, ei, eoo)'s (modulo Z) we should try is given in [33, 
Table (I), Section 4]. To make the expression size as small as possible, 
the (eo, 61,600) that leads to the lowest degree for / should be tried 
first. The table shows how the degrees are related to each other; to 
minimise the degree, we should check (0,0,0) first, and (0,1/2,1/3) 
last. The first case (0, 0, 0) did not lead to a solution. The second case 
we tried was (0,0,1/3), and since that led to a solution, we did not 
need to try more cases. The table in [33] shows some of the other cases 
must also lead to a solution; one sees in the table that (0, 0, 1/3) (and 
hence P2) is solvable in terms of (0, 1/2, 1/6) and (0, 1/2, 1/3). We can 
also read from the table that if we had used those exponent differences 
for our 2F1, then the degree of / would have been 2, respectively 4, 
times higher. 

A Maple session to perform all operations of this section is given 
online at http://algo.inria.fr/chyzak/Rooks/Compute2Fl.mpl. 

3. Further comments 

3.1. Asymptotics. It was pointed out in [16] that general results 
from [27] on the asymptotics of diagonal coefficients of multivariate 
generating functions imply 

, , 9^3 64" 

407r n 

Here we show that the same result can be derived in a simpler way 
(using univariate, instead of multivariate, singularity analysis), based 
on the explicit expression (16). 

Indeed, basic singularity analysis [19, 20] shows that, since the dom- 
inant singularity of dx(Dia,g{f)) is x = 1/64, with residue 

ii2-m _ 32 ^ /1/3 2/3 



6 /1/3 2/3 

4T'2^M 9 



64^ 




5 ^^'1 
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the asymptotic behaviour of the coefficient hn of dx(Di8ig{f)) is hn 

r 64"^ 

r ■ 64", and thus a„ ~ — ■ . 

64 n 

The value of r is found using Gauss's formula 



2F1 



a b 



T{c)T{c-a-b) 
T{c-a)T{c-b) 



and the multiplication formula 



riz)-r{z + -]---r{z 



m 



m 



with m = 3 and z = 1. This gives r = 72^/3 / (dn) , and proves (18). 

3.2. Non-algebraicity and simpler formulas. One may wonder 
whether the function Diag(/) is algebraic. This is not the case, and 
can be seen in several ways: (^z^ as a consequence of the asymptotic be- 
haviour of the coefficients a„ (generating series of sequences asymptoti- 
cally equivalent to n"p" with a G Zi_ are necessarily transcendental [18, 
Thm D]); (ii) by doing a local analysis of P near its singularities (the 
presence of logarithmic terms in a local expansion excludes algebraic- 
ity); (in) by applying Schwarz's classification [28] of algebraic 2-^1's. 

It is also legitimate to seek for an even simpler closed form, not 
involving any integration sign. Such a reduction is not possible: if it 
were, the algorithm given in [2] should have found it. The sequence (a„) 
satisfies the third-order recurrence (15), but no second-order recurrence 
with polynomial coefficients, for otherwise (15) or a closely related 
adjoint recurrence would need to possess a hypergeometric solution, 
which is easily sorted out by Petkovsek's algorithm [26]. 

3.3. An alternative hypergeometric expression. Using a com- 
pletely different method. Frits Beukers [8] independently found that 
the diagonal 3D rook path series G{x) satisfies 

1 



(19) 
where 

with 



G'i 



X 



x] 



2{l + 6x) 



Ax)H'[ 



X] 



4H x) 



Hix) 



1/4 

92 



'1/12 5/12 



92= {I- 4x)(l - QOx + 120x^ - 64x^), 

J{x) 



J{x) 



9l 



1728(1 - xfx^{2 - 3x)3(l - 64x) ' 
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The function J{x) is obtained, by a method similar to the one in [7], 
as the J-invariant of a family of elliptic curves naturally associated to 
the integral interpretation of the diagonal problem. 

Of course, once obtained, it is an easy task to check that the new hy- 
pergeometric expression coincides with ours: the difference of the two 
hypergeometric expressions satisfies a simple linear differential equation 
derived from the defining linear differential equation for the hypergeo- 
metric 2-^1, and it is easily expanded as a Taylor series, showing that 
enough first terms are zero, so that the difference is zero. This calcula- 
tion is available as Maple code at http : // algo . inria . f r/ chyzak/Rooks/ CompareFormulas , 

Interestingly, one can algebraically reduce one expression to the other 
by combining only two classical hypergeometric identities: 

(a) Gauss' contiguity relation 



(20) 



9(1 



X) 



■ 2 



'1/3 2/3 



obtained by evaluating at a 
general contiguous transformation 



'1/3 2/3 



|, and c 



X 



1 the more 



c(l 



X 



d a h 



X 



dx 

c-a){c-b) ■ 2F1 



a b 
c+ 1 



x] + c{a + b - c) ■ 2F1 



(b) Goursat's quartic relation 



(21) 
2F1 



'1/3 2/3 



X 



iFi 



'1/12 5/12 



64x^(1 



x] 



(9 



ix 



a b 
c 



X 



-1/4 



obtained by taking a = ^ in the identity [21, Eq. (126)] (see 
also [37, Eq. (25)]): 



iF, 



, 6a + i 



X 



iF^ 



a a + I 
2a + I 



QAx^{l-xY 
(9 - 8x)3 



^ ^x 
~"9", 



3a 



The detailed reduction is as follows: Equations (20) and (21) show 
that 



iFi 



'1/3 2/3 



9(1-^ _ / ^ /1/12 5/12 



18 A. BOSTAN, F. CHYZAK, M. VAN HOEIJ, AND L. PECH 

where ip{x) = 64x^(1 — a:)/(9 — Sx)^ and xi^) = ~ ^) ^ • By 
substitution, the latter equahty imphes 

'1/3 2/3 



2 



ip{x) 



9 1 -^(x) / ^ /1/I2 5/12 



with V5(x) = 27x(2 — 3x)/(l — 4x)^. The right-hand side term rewrites 



;i -64x)(l -x)(l -4x) / ^ A/12 5/12 



12(1 + 6x) V V ^ 



1 \ / \ -V4\ ' 



j{x)j \{i-Axy 



and this, in combination with equation (16), finally implies equation (19). 

3.4. The case of 3D queens. A closely related problem is to find the 
generating function for the diagonal sequence for 3D queen paths. This 
is Sequence A143734 in Sloane's encyclopedia http : //oeis . org/A143734, 
whose first terms read 

1, 13, 638, 41476, 3015296, 232878412, 18691183682, 1540840801552, . . . 

Computationally, it turns out that this problem is more difficult 
than the one on 3D rooks. By Lemma 1, it boils down to finding the 
complete diagonal Q{x) of the rational function 

S t X st tx xs stx 



1 — s 1 — t 1 — X 1 — st 1 — tx 1 — xs 1 — stx 

We empirically found (using guessing routines based on Hermite- 
Pade approximants) that Q{x) satisfies a differential equation Lq{x, d^) 
of order 6 and polynomial coefficients (in x) of degree at most 71. The 
corresponding sequence satisfies a recurrence of order 14 with polyno- 
mial coefficients (in n) of degree at most 52. Moreover, it is very likely 
that these are the minimal-order equations. They are available online 
at http://algo.inria.fr/chyzak/Rooks/QueensRecEq.mpl for the 
recurrence and http : //algo . inria.fr/chyzak/Rooks/QueensDiffEq.mpl 
for the differential equation. 

The conjectured equations already have a large size, so that proving 
their correctness by the creative-telescoping approach becomes a chal- 
lenge; we have not yet succeeded in proving this in a fully algorithmic 
fashion. 

However, we are convinced that our guessed equations are correct: 
they pass all the various check tests described in [10]. In addition, the 
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degree-71 leading coefficient of the differential equation factors as 

■ (5x - 4) ■ (a; + 1) ■ (a; - 1)^ ■ (SGx^ - 40x + 9) ■ 
(512x^ - mix^ + 84x^ + 95x - 1) ■ 
(14063x^ - 15940x^ - 5918x^ + 12063x^ - 4118a;^ + 575x - 40) ■ P(x) 

where P{x) G Q[x] is irreducible of degree 49. 

Note the known fact [27, 16] that the nth coefficient of Q{x) is asymp- 
totically equal to Kc~^"/n, for some k G M and c ~ 0.2185. More 
specifically, c is the smallest positive real root of the equation 

\l)x-l \2)l-x^ \?,)l-x^ 

which is consistent with the presence of the factor 512x'^ — QQlx^ + 
84a:^ + 95x — 1 having ~ 0.0104 as a root (in fact, its root that is 
closest to zero). 

The sixth-order differential operator LQ[x,dx) can be factored as 
Lq = L^- r['^^ ■ R^^ , where L4 has order 4, and both R^""^ and R^^ have 
order 1. By analogy with our result on rook paths, a natural question 
is whether L4 is solvable in terms of indefinite integrals, algebraic func- 
tions, and 2F1S. expressions. If it were, then it would be solvable in 
terms of solutions of second-order equations, and by [29] at least one 
of the following should then be true (see [35, Sections 1.1 and 1.2] for 
definitions and details): 

(1) L4 is reducible, 

(2) L4 is gauge equivalent (under Maple's DEtools[Homomorphisms]) 
to the third symmetric power of a second-order operator, 

(3) L4 is gauge equivalent to the symmetric product of two second- 
order operators. 

Regarding option (1), the operator does not factor (D Factor in Maple). 
Regarding options (2) and (3), [35, Section 3] shows how to find order 
reductions of such type (if they exist). We found that no such reduction 
exists. Option (2) can be ruled out quickly by local data (if a third 
symmetric power has a logarithmic solution at x = p, then the cube of 
a logarithm should appear as well). For option (3), we computed the 
eigenring of the second exterior power of L4. To conclude, the fourth- 
order factor L4 is not solvable (in the sense of Singer's definition [29]) in 
terms of solutions of second-order equations; in particular, this means 
that, in this definition, there is no expression in terms of 2-^1 functions 
for the generating function Q{x) of diagonal 3D queen paths. 
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It appears that L4 has no 4F3 type solution either, but we have not 
obtained a proof of this fact yet. 

3.5. Faster but heuristic approach to getting the certificates. 

As already mentioned in the introduction, we have chosen in Section 2 
to confine ourselves to calculations that follow an algorithmic approach, 
that is, that would provably (theoretically) succeed on any other ra- 
tional input. Beside this, heuristics are available to cleverly guess a 
solution to the key equation, including certificates, and often lead to 
faster computations when they succeed, although they can fail in the- 
ory and sometimes do in practice. One of the most advanced strategies 
in this spirit is due to Koutschan [23]. Note that by rational- function 
normalisation, verifying the key equation a posteriori proves the ob- 
tained results. 

We have implemented such a strategy in the rational case, and we 
could get a solution (P, S",T'), where P is the same operator as in 
Section 2.4.1, but for other certificates S' and T'. In effect, this heuris- 
tically solves a partial differential equation for (some of) its rational 
solutions. The calculation is 4 times faster than in the algorithmic way. 
Interestingly enough, the denominators of S' and T' only involve qi, 
and not its discriminant disct(gi). 

Acknowledgements. We wish to thank the referee for his quick and 
serious job and for his constructive comments on the manuscript. 
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Appendix A. Details of calculations of the iterated 

CREATIVE TELESCOPING FOR THE 3D ROOK 
GENERATING SERIES 

In this appendix, the calculation sketched in §2.4.2 is explained in 
more detail, and we provide explicit values of intermediate computa- 
tions. We start by introducing the notions needed for our description. 
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Differential operators. Whereas the sketch of the algorithm was given 
with references to the function F to be integrated and its derivatives, 
here we prefer presenting the calculations as they are done in actual 
implementations: a primal role is played by differential operators that 
represent linear combinations of derivatives of a given function. 

As certain operators map the function under consideration to zero, 
operators have to be viewed up to addition of such annihilating oper- 
ators. This explains the use of congruences in the calculations below, 
denoted Li = L2 mod /, when Li{F) = L2{F) and / is the (left) ideal 
of all operators that annihilate F. 

Representing a D-finite function. Recall from [30, 25] the definition 
that a D-finite function of variables u and w is a function F for which 
there exists a finite maximal set of derivatives d^d^{F) that are linearly 
independent over Q(m, f). The natural question is how a function is 
algorithmically recognised and represented as D-finite. 

The best algorithmic answer would be in terms of non-commutative 
Grobner bases (and ideals of Hilbert dimension 0) [12], but here, to keep 
the sophistication to a minimum, we shall use the alternate definition 
that a function is D-finite if there are non-zero operators A{u,v,du) 
and B{u,v,dy) that annihilate it, one for each derivative. Follow- 
ing [12], the pair {A, B) of annihilating operators will be called a rect- 
angular system. 

Rational solutions of parametrised systems of linear differential equa- 
tions. Chyzak's algorithm relies on an algorithm that, given matrices 
A G X„(K(t)) and B G MnxdiJ^if)) for some field IK that depends on 
the application, solves 

(22) ^ + Ay = Be 

for y G K(t)"' and e G W^. In our uses, the field K is a rational 
function field in the other variables of the problem. We do not describe 
this algorithm here, and refer the reader for instance to [1] (and the 
description of its use in [11]). 

For convenience, let us recall the notation of §2.4.1. Let q be the 
denominator of F, then q = stqi with qi = —st + 2s^t + 2t^ + 2xs — 
3st^ — 3xt — 3xs^ + 4xst an irreducible polynomial of Q[x, s,t]. Let 
disct(gi) = (x — s)(16xs^ — 4s^ — 24a;s-|-4s^-|-9x — s) be the discriminant 
of qi with respect to t. 
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A.l. Stage A: First iteration of Chyzak's algorithm. Stage A 
in §2.4.2 looks for operators P^"^ G Q{x, s){dx,ds) (according to the 
ansatz (10)) and rational functions 0^"-* G Q{x, s, t) such that P^"\F) = 
f (0(°)F) (cf. (11)). A property of the left ideal (P(°)) thus obtained 
in Q(a;, s){dx, ds) is that it contains a rectangular system. 

Ansatz with r = 1 . We first look for operators P'-"^ of total degree r = 
1. Let ?7o,05 ■'71,0, ''70,1 £ Q{x,s) be unknowns such that P = ?7o,o + 
^i,o^rr+'7o,if^s5 and let G Q{x, s, t) be another unknown. The equation 
p(F) = f (0F) rewrites 

dF dF d(t) dF 

which is an equation in and t]q^q, ?7o.i, ''7i,0; of the form (22) with K = 

Q(x, s). All solutions are proportional to the one given by 

(23) 

r/o,o = 2(s-l)(3s2-6s + 2), 

= Qxs^ - 2s^ - lOxs^ + - Ax'^s + Wxs + 3a;^ - Ax, 
r/o,i = 2s(3s-2)(s-l)2, 

^ _ -t{s - l){-6xs^ + QsH - + llxs - 9st - Ax - xt + At) 

t — X 

Set Pi to the operator P obtained for these specific values. The ideal 
generated by it obviously does not contain any rectangular system in 
Q{x,s){d,,d,). 

Ansatz with r = 2. Since all the monomials dl can be reduced by Pi, we 
set r = 2 and look for P = rjQ+rjidx+rj2d1, with ?7o, rji, 772 G Q.{x, s), and 
4> G Q(x, s,t) such that P{F) = ^(0F). (This ansatz is not strictly 
of the form (10) that would involve ds] this will require additional 
manipulations when obtaining a rectangular system.) The equation 
on F becomes 

OF d^F d(j)^ dF 

VoF + VijT + = 7^^ + 

ox ox^ at at 

which is a parametrised equation in and 770, Tii,ri2, with K = Q(x, s). 

All solutions are proportional to the one given by: 

^70 = 0, 

r/i = -2(-19s^ - 9a; + 13s^ + 7s - IQxs^ + 2Axs), 
(24) \ ri2 = discf(gi), 

^ -t{3s - 2)(s - l)(2s2 - 4st - s + 3t)(s - t)^ 

*^ = 77 ^ ' 

[t - x)qi 



EXPLICIT FORMULA FOR DIAGONAL 3D ROOK PATHS 



25 



and we call P2 the operator P obtained for these specific values. The 
ideal {Pi, P2) now contains a rectangular system: eliminating from 
Pi and P2 results in an operator = Vo ~^ v'l^s + V-^^s ^ with coefficients 
of respective degrees (4,9), (4, 10), and (4, 11) in (x, s). This ends the 
calculation. 

The Maple code to compute P2 can be obtained at http : / / algo . inria . f r/ chyzeik/Rooks/Re 

A. 2. Stage B: Second iteration of Chyzak's algorithm. To per- 
form Stage B in §2.4.2, we now work modulo the ideal / = (Pi, P2) of 
Q(x, s){dx, dg). The vector space V = Q(x, s){dx, dg) / 1 is 2-dimensional 
over Q{x, s), with basis (1 + I,dx + I). In other words, we consider a 
generic solution F of Pi and P2, and represent it by 1 + /. 

To simplify the notation in what follows, we normalise the operators 
into the form Pi/(2s(3s — 2)(s — l)^) = ds—pidx—po and P2/disct(gi) = 
dl — qdx, for coefficients po,Pi, q G Q(x, s). 

The algorithm now starts from ansatz (8) for P, but with the form (9) 
for S, in order to find a relation of the form (12). For increasing values 
of d, it looks for an operator P = X^iLo Vi^^x^ such that P = dgQ mod /, 
with rii G Q(x), 0o, 0i e Q(x, s), and Q = 0o + (t^idx- 

There is no solution for (i < 3, and we give details for c? = 3 only. 
First, the algorithm simplifies (7). Modulo J, we have = qd^, and 
therefore 

Additionally, dg = po + pid^, so 

OX \ OX J Ox \ Ox J 

As to the right-hand side: 

dsQ = ^ + (j)ods + + (f)idsd^ 

OS OS 

, f , ^ \ A , ^^^^ , A (^Po , ( , , ^ ;^ ^ 

= ^— + [po + pic^x) 00 + -K-Ox + \ P0 + ^+Piq]Ox] ■ 

OS OS \ Ox \ Ox J J 

Since P = r]o + (r]i + qri2 + (|f + 9^) ^73) ^^x, the congruence P = 
mod / expressed in the basis (1 + I,dx + I) is equivalent to the differ- 
ential system 
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ds 



ds 



-Po(Po - ^(Pi + Vo, 
ax 



-Pi4>o 



Po 



dpi 
dx 



Piq 01 + 77i + (1V2 




Solving this parametrised system of hnear differential equations in \l{x) 
results in a 1-dimensional vector space, with basis given by the solution 



(26) 



f ^0 

V2 
V3 



0, 

4(576x2 



801x' 

„2 



108x + 74), 
4608x^ + 813x2 - 6372x2 + 514x - 4, 
x(x - l)(64x - l)(3x - 2)(6x + 1), 
(53 + 108x)(3s -2)s 

7 



2(s - l)MisCi(gO 

where 7 is an irreducible polynomial in Q[x, s] of degree (5, 7) in (x, s). 

The Maple code to compute 7 can be obtained at http : // algo . inria . f r/ chyzak/Rooks/Sec 

A. 3. Stage C: Reconstruction. We finally derive an explicit identity 
of the form (13) in Stage C in §2.4.2 and obtain the corresponding 
values for S and T given by (14), thus leading to the solution (P, S, T) 
of (3). 

The previous calculation of Pi and P2 can be reinterpreted as follows: 
Introducing the left ideal J that annihilates F, that is, the ideal J 
generated by 



Fa -'^l 

ax 



Fd-"^ 

as 



dt ' 



we have obtained the congruences 



Pi = dtipi mod J and P2 = dt4'2 mod J, 

for rational functions tpi and ?/'2 th.a,t are the functions in (23) and (24), 
respectively. Additionally, P = dgQ mod / for P = 770 + rjid^ + 772^^ 
and Q = 00 + 4>idx described by (26), and I = (Pi, P2). 

To express P — dsQ G (Pi, P2) more explicitly, we look for Ai,A2 G 
Q{x,s){dx,ds) such that P - dsQ = AiPi + A2P2. This division is 
made algorithmic by observing that (Pi,P2) is a Grobner basis, and 
by performing a reduction by this Grobner basis. In simple terms, the 
calculation is as follows: We first use Pi to eliminate dg from P — dsQ, 
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then use P2 to divide the remainder. Collecting the quotients, we then 
find: 

108X + 53 71 

' 2(s-l)3 4s(3s-2)(s-l)Misct(gi) 

^ = 72 

^ 4s(3s-2)(s-l)Misc4(gi)2 

x{x- l)(64x- l)(3x-2)(6x+ 1)^ 

disct(gi) 

with 7i, 72 irreducible polynomials in Q[a;, s] of respective degrees (5, 7) 
and (7, 10) in (x, s). 

Since -Pi(-F) = dt{ipiF) and -P2(-f') = dt{ip2F), we get the equality 
(P-(9,Q)(F) = (v4i(9tVi + ^9tV'2)(i^)- Since and A2 do not depend 
on t, we get P(F) = {dsQ + dt{Aiij, + A2^2)){F). 

Let S = Q{F) and T = (Ai^^i + v42^/'2)(P). Then S,T e Q{x,s,t) 

dS _i_ dT 
ds dt 

form 



and -P(-F) = ^ + K turns out that S and T are of the very special 



g_ {s-t)U ^_ {s-t)V 



with U,V irreducible polynomials in Q[x, s,t] of respective degrees 
(5,8,3) and (8,14,5) in {x,s,t). 

The Maple code for this reconstruction of S and T can be obtained 
at http : // algo . inria . f r/ chyzak/Rooks/Reconstruct ion . mpl . 
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